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Controlling semi-convergence 
phenomenon in non-stationary 
simultaneous iterative methods 
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Abstract 


When applying the non-stationary simultaneous iterative methods for 
solving an ill-posed set of linear equations, the error usually initially de- 
creases but after some iterations, depending on the amount of noise in the 
data, and the degree of ill-posedness, it starts to increase. This phenomenon 
is called semi-convergence. We study the semi-convergence behavior of the 
non-stationary simultaneous iterative methods and obtain an upper bound 
for data error (noise error). Based on this bound, we propose new ways 
to specify the relaxation parameters to control the semi-convergence. The 
performance of our strategies is shown by examples taken from tomographic 
imaging. 


Keywords: Simultaneous iterative methods; Semi-convergence; Relaxation 
parameters; Tomographic imaging. 


1 Introduction 


A mark-point in the history of medical imaging, was the discovery by Wil- 
helm Rontgen in 1895 of x-rays [10,22]. The problem of generating medical 
images from measurements of the radiation around the body of a patient 
was considered much later. Hounsfield patented the first CT-scanner in 1972 
(and was awarded, together with Cormack, in 1979 the Nobel Prize for this 
invention). This reconstruction problem belongs to the class of inverse prob- 
lems, which are characterized by the fact that the information of interest is 
not directly available for measurements. The imaging device (the camera) 
provides measurements of a transformation of this information. In practice, 
these measurements are both imperfect (sampling) and inexact (noise). 
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Figure 1: Semi-convergence phenomenon 


The mathematical basis for tomographic imaging was laid down by Jo- 
hann Radon already in 1917 [20]. The word tomography means “reconstruc- 
tion from slices”. It is applied in Computerized (Computed) Tomography 
(CT) to obtain cross-sectional images of patients. Fundamentally, tomo- 
graphic imaging deals with reconstructing an image from its projections. The 
relationship between the unknown distribution (or object) and the physical 
quantity which can be measured (the projections) is referred to as the forward 
problem. For several imaging techniques, such as CT, the simplest model for 
the forward problem involves using the Radon transform R, see [1, 16, 18]. 
If x denotes the unknown distribution and 8 the quantity measured by the 
imaging device, we have 

Rx = 6. 


The discrete problem, which is based on expanding y in a finite series of 
basis-functions, can be written as 


Ax ~ b, (1) 


where the vector b is a sampled version of 6 and the vector x, in the case 
of pixel-(2D) or voxel-(3D) basis, is a finite representation of the unknown 
object. The matrix A € R™*", typically large and sparse, is a discretization 
of the Radon transform. An approximative solution to this linear system 
could be computed by iterative methods, which only require matrix-vector 
products and hence do not alter the structure of A. 

Initially the iteration vectors approach a regularized solution while con- 
tinuing the iteration often leads to iteration vectors corrupted by noise. This 
phenomenon is called semi-convergence by Natterer [18]; for analysis of the 
phenomenon, see, e.g., [1,2,9,11,13,19,21]. The typical overall error behavior 
is shown in Figure 1. 

The Algebraic Reconstruction Technique (ART) is a fully sequential 
method, and has a long history and rich literature. Originally it was pro- 
posed by Kaczmarz [15], and independently, for use in image reconstruction 
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by [13]. The vector of unknowns is up-dated at each equation of the system, 
after which the next equation is addressed. In the simultaneous algorithms 
the current iterate is first projected on all sets to obtain intermediate points, 
and then the next iterate is made by an averaging process, as convex com- 
bination, of intermediate points. The prototype of these algorithms is the 
well-known Cimmino method [5]. We now explain block-iterative method. 
The basic idea of a block-iterative algorithm is to partition the data A and 
b of the system (1) into blocks of equations (rows) and treat each block ac- 
cording to the rule used in the simultaneous algorithm for the whole system, 
passing, e.g., cyclically over all the blocks, see Figure 2. 

An iteration vector of the non-stationary simultaneous iterative method 
(SIM) is defined as follows 


Lk41 = tE+A,ATM(b— Azz), &=0,1,--- (2) 


with zo € R” where {A,}?2, are relaxation parameters and M is a given 
symmetric positive definite (SPD) matrix which depends on the particular 
method. In some papers in image reconstruction from projections, the term 
“simultaneous iterative reconstruction technique (SIRT)” is used for “SIM”; 
see, e.g., [7,8, 21]. Several well-known simultaneous methods can be written 
as (2) for appropriate choices of the matrix M. With M = I we get the 
classical Landweber method [17]. Choosing M = +diag(1/||a;||?) where a; 
denotes the ith row of A leads to Cimmino’s method [5]. The CAV method [4] 
uses M = diag(1/ 1 Nja;;) where Nj; is the number of non-zeroes in the 
jth column of A. 

We study semi-convergence behavior of the non-stationary SIM, when 
applied to noisy data. Our main focus is to propose some techniques for 
updating relaxation parameters to control the data error. Having a reliable 
stopping rule leads to stop the iterative method in an iteration which makes 
a proper approximation of the sought solution. Otherwise, we may stop 
the iteration process early or far from a proper iteration index. For this 
reason we introduce relaxation parameters to postpone the semi-convergence 
phenomenon. 

The iteration index of an iterative method may be considered as a regu- 
larizing parameter. We explain this a bit more. Let x* be the sought solution 
using exact data and let £, and x, denote the iterate using noisy and exact 
data respectively. Then we have 


Ze — 2*|| < [Ze — eel] + [lax — w" I]. (3) 


Therefore, the error decomposes into two components, the data error (or noise 
error) and the approximation error (or iteration error). The semi-convergence 
of the iteration interplays between these two error terms. 

The semi-convergence behavior of the SIM with constant relaxation pa- 
rameter is analyzed in [8] where the related M-matrix is a symmetric positive 
definite (SPD) matrix. Based on this stationary, they suggest two strate- 
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Figure 2: (right to left) sequential method, simultaneous method and sequential block- 
iterative method 


gies for picking relaxation parameters to control the upper bound of data 
error. The obtained sequence of relaxation parameters is nonnegative and 
nonascending. 

Later in [7], the projected version of the non-stationary SIM is studied 
where the M-matrix is again assumed SPD. As [8], they consider nonascend- 
ing sequence of relaxation parameters and emphasize both strategies of [8]. 
In [7], using nonexpansivity of the projection operator leads to assuming two 
cases, i.e., the full column-rank problem (rank(A) = n) and the rank-deficient 
problem (rank(A) < n) which is handled by a slightly modified problem. 
Furthermore, they present upper bounds for noise error and iteration error 
where rank(A) =n, see [7, Theorems 3.3 and 3.8] respectively. Also those 
bounds can be achieved for the modified problem with an unknown regular- 
ization parameter (see [7, (3.22),(3.23)]) under some assumptions [7, Lemma 
3.9]. In section 2, we give an analysis of non-stationary SIM without having 
any restriction on rank(A). Additional to strategies given in [8] and [7], we 
introduce another strategy for choosing relaxation parameters which is able 
to make more reduction in noise error upper bound comparing with the old 
strategies. 

In Section 3, we consider SIM and give its semi-convergence analysis with 
three strategies for picking relaxation parameters. We demonstrate the per- 
formance of our strategies by examples taken from tomographic imaging in 
Section 4. 


2 Simultaneous iterative algorithm 


In this section we give an analysis of the non-stationary SIM without assum- 
ing any restriction on rank(A). 

Let ||a|| = Va? x and |la|| 4. = Vz? Ma denote the 2-norm and a weighted 
Euclidean norm respectively. Also, let M!/? and p(Q) denote the square 
root of M and the spectral radius of Q respectively. For W € R™*”, we 
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use N(W) and R(W) to denote the null space and range of W respectively. 
The orthogonal projection from R” onto N(W) is denoted by P(W). Also 
the orthogonal complement of a subspace K of R” is denoted by K+. Here 
«tm (A,b) denotes a solution of min || Ax — b||;¢ with the minimal Euclidean 
norm. 


The convergence analysis of SIM can be obtained in, e.g., [14, Theorem 
II.3] and [3]. 


Theorem 1. Let p = p(A?MA) and assume that 0 < € < Ax < (2—6)/p. 
Ife >0, ore =0 and XP, min(pAx, 2 — pAx) = 00, then the iterates of (2) 
converge to x(A,b) + P(A)xo0. 


2.1 The error in the k-th iteration 


As we mentioned before, in this section, we give the same upper bound as [7, 
Theorems 3.3 and 3.5] but without any restriction on rank(A). Based on 
our analysis, we give another strategy for choosing relaxation parameters. 
This strategy is capable to reduce noise error upper bound more than the old 
strategies given in [8] and [7]. 


Let B = A?MA , x* = 2,(A,b) and consider the singular value decom- 
position (SVD) of M'/?A as 


M124 =USVT 


where © = diag(o1, ...., %p,0,...,0) € R™*” with a) > 02 >... > 0, > 0 and 
p is the rank of A. Let z, = x, — a*. Using (2) we have 


Zk41 = Zk + ATM (b — Az, - Ax”) = (I = AnB) Zr 


which leads to 
k-1 


Zk [[v _ Ap—1-iB)2Z. 


i=0 
Since zg) = Xo — 2”, we obtain 
0 0 ’ 


k-1 


ay =a" + |] (2 An-1-1B) (a0 - 2”). (4) 
1=0 


Using the orthogonal decomposition theorem, we have R” = N(B)@ N(B)- 
and N(B)+ = R(B). Therefore we get 29 = £9 + P(B)xo where %o € N(B)+ 
and P(B)xo € N(B). Thus we can rewrite (4) as 
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k-1 
ap =a" + [] (2 — An-1-iB) (@o + P(B) a — 2"). (5) 
i=0 
Since BP(B)ap = 0 and P(B) = P(A), we obtain 
k-1 
[[U = Ax-1-:B)P(B)20 = P(B)x0 = P(A)20. (6) 
i=0 


Since #9, x* € N(B)+, we can rewrite # — x* as 
P 
Lo -xf= ; CjU5 (7) 
j=l 


where c; and v; are scalar and the j-th column of V respectively. Using (5), 
(6), (7) and the SVD of B, we obtain the following expression for the k-th 
iteration 


k-1 p 


ap = a* + P(A)xo + | [ (2 - Ana-iVETZV™)(S cjvj) = 2* + P(A)z0+ 
i=0 j=l 
k-1 k-1 p 
+ Vdiag (Tle = dio7), tansy [[a = Xio3), 1, sans : VTS | ¢9; 
i=0 i=0 j=l 
p k-l 
= 2" + P(A)ao t+ ¥> [[ 1 —rio)ejay. (8) 
j=1 i=0 
Let 6 = b+ 6b and 
Trt = Te+ dp ATM (b _ AZ;) (9) 


where 60 is the perturbation consisting of additive noise. Setting 7, = ©,—2"*, 
we get 


Zk = Ze-1 + dn—-1 47 M(b + 6b— AZp_-1 — Ax*) 
= (I — \p_1 BYZp_1 + Xg_-1 A? Mb 


k-1 k-2 k-1 
= [[v — Xp-1-iB)20 + S- II (I — An—-1-jB)\ AT Mbb+ 
i=0 i=0 j=itl 
+ Ap—1 AT Mob. (10) 


Since %p = Xo, similar to (8), we have 
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p k-1 
tp= 2 +P Ajto+ S- [[a - NiO} ejuj+ 
j=1 i=0 
k-2 k-1 
+ 5° [J = An-1-5B)AAT Mb + Ap_1 AT Mb. (11) 
i=0 j=it1 


We now assume that the sequence of relaxation parameters is nonnegative 


and nonascending, i.e., 
0< Au SA: (12) 


and consider the following function introduced in [8] 


1— (1 —20?)* 
ee 


U*(o,A) = (13) 


Theorem 2. Let w = ||M1/?6b|| and 0 < Ax < +e. The noise error of SIM 
ak 
is bounded above by 


wrAQo1 


|Z. — cell < Toray ep Ae-1): (14) 
— Pp 


Proof. By subtracting (8) and (11), we obtain 


k—-2 k-1 
Ee-te= >) |] (Ana B)AiA bb + Ag_1A™ MOD. (15) 
i=0 j=itl1 


Therefore we have 


k-2 k-1 
Ze —aell < So AG| [] (—Ax—a-7B)A™ M60]| + Ax—a||A7 MOD]. (16) 
i=0 j=it+l 


Using the SVD of M'/2.A, we get that 


k-1 k-1 
II (I — Ax-1-;B)AT M1? =V II (I —Apa—phTU)VIVET UT 
j=itl j=itl 
=VW;,.U7 
where 
k-1 k-1 


Win =diag{ [J] -Ajof)or,..., [] (-Ajo?)op,0,...,0 
j=itl j=itl 
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Using (12) and 0 < Ag < sy, we obtain 
1 


k-1 k-1 
my : T ygi/2 se ve mena 
TD f= Xe-1-5B)A7 M7 < [el] = max | T] (1 -As03)0s 
j=itl g=itl 
< o1(1 = Aedes) (17) 


Since ||A?M6b|| < o1w, we conclude that, using (12),(17) and the assump- 
tions of theorem, 


k-2 k-1 
|Z 2. = Lp|| < S- Agw II (I = Ap—1-7B)AT MV? + Ago iw 
1=0 j=it+l 
2 


< oS Anwar (1 = Nii) + Ago iw 


s=0 
WwArAQo1 1- (1 . Ar—102)* 


Ak-19p Op 


WAQo1 & 
= ———_ V Ak-1): 
An-19 p (op, k 1) 


This completes the proof. 


Remark 1. To obtain a similar result as (14) where the projected case of (2) 
is employed, we refer to [7, Theorem 3.3] where it is assumed rank(A) =n. 
Similar to [8], we consider the equation 


ge—1(y) = (2k — Ly*-* — (yk? +--- ty +1) =0 (18) 


which has a unique real root ¢, € (0,1). The roots satisfy 0 < Cy < Ck41 <1 
and limg_-.o Ck = 1 (see [8, Propositions 2.3, 2.4]), and they can easily be 
precalculated, see Table 1. 


Again, let a; denote the largest singular value of M‘/?.4. Then we have 
the following alternative upper bound for the noise error. 


Theorem 3. Assume that a1 < 1/./XAp-1; then 


WAgC1 1- a 


V Ak—-198n V i Ce 


where Cy, is the unique root in (0,1) of (18). 


Ilvx — Fell S (19) 
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Table 1: The unique root ¢, € (0,1) of gx—1(y) = 0, ef. (18), as function of the iteration 
index k 


0.8392 | 13 0.9019 | 18 0.9294 | 23 0.9449 | 28 0.9548 
0.8574 | 14 0.9090 | 19 0.9332 | 24 0.9472 | 29 0.9564 
O 0.8719 | 15 0.9151 | 20 0.9366 | 25 0.9493 | 30 0.9578 
1 0.8837 | 16 0.9205 | 21 0.9396 | 26 0.9513 | 31 0.9592 


k 

0.3333 7 0.8156 | 12 0.8936 | 17 0.9252 | 22 0.9424 | 27 0.9531 
8 
9 


aookwn sr 
So 
ol 
ol 
io) 
w 


Proof. Using [8, Proposition 2.3] we obtain the following bound for the func- 
tion W*(c, \) appearing in (14): 


max WU (ai, Ag—1) < max U* (a, Ak—1) 


l<i<n O0<o<o, 
< max WU (a, n-1) < ‘(a mas (20) 
~ O<a<1/4/Anaa ~ V1—¢k 


The assumption in the theorem implies 


oy <1/VrAr-1 S Ni a (21) 
Then by (14) and (20), and assuming (21), we obtain the bound in (19). 


Remark 2. The case Ax—1 € (1/07, 2/07) is discussed in [8, Remark 2.2]. 


3 Choice of relaxation parameters 


Using (19), we propose following strategies for choosing relaxation parame- 
ters: 


42 for k =0,1 
Wi, -—rule: Ap= i (22) 
2 
sr(1—<¢,), fork > 2, 
al 
ney fork =0,1 
Wo—-rule: A= : (23) 


2(1-G)(-¢h)-?, for k> 2 


W3-rule: Ap= (24) 


ee for k = 0,1 
1 
Ze(1— Ce)’ 1(1 — Cf), for k > 2 


where {¢x}x>2 are the roots of (18) andl <r< 2. 
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Noise error upper bound 
T T 


Figure 3: Noise error upper bound (25) for different strategies with the factor w/op 
omitted 


Remark 3. Using (19) and strategies (22-24), we have the following upper 
bounds for noise error 


ee eG) en 
[Ze — eell < ¢ SA Cf)P(L— Ce), We (25) 
a a, Gy, V3 


for k > 2. 


Figure 3 shows the behavior of noise error upper bound (25) for different 
strategies. As it seen, V3 with r = 1 and W3 with r = 3 give the smallest 
and largest upper bounds respectively. Furthermore, V3 with r = 1.5 gives 
smaller upper bound than WV; and Wo. 


Remark 4. It is easy to check that, using [8, Theorems 3.1 and 3.3], the 
both strategies (22) and (23) satisfy all conditions of Theorem 1. Therefore, 
the sequence x, generated by (2) converges to x(A, b) + P(A). 


Next we will check that the relaxation parameters defined in (24) satisfy 
all conditions of Theorem 1. 


Proposition 1 The sequence generated by (2) with strategies (24) converges 
to xy (A, b). 


Proof. Since p = 07, we have 0 < pAx < 2. Using [8, (2.17), (3.10)], we obtain 
that 
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ae ie k-1\7 
1 r—lyy _ ¢ky2 joe ee 5 eee 
SAG ate Ik+1 1 
k>2 k>2 
1 ra k 2 
-¥ (sn) (so) 


k>2 


k2 
ee (26) 
a (2k + 1)r+ 


It is clear that (26) diverges if r < 2. Therefore, we have Ux>2Ax~ = 00. It is 
easy to check that min{pA,z,2— pA} = pAx for k sufficiently large. Thus, all 
conditions of Theorem 1 hold and consequently the sequence x, generated 
by (2) converges to xy¢(A, b) + P(A)zo. 


4 Numerical results 


In this section we give two examples of computerized tomography field. We 
used 5% and 10% white Gaussian noises to produce noisy data. The constant 
optimal relaxation parameter Ajp¢ refers to the strategy when a constant 
value of the relaxation parameter is used, chosen such that it gives rise to the 
smallest relative error within 20 iterations. For the choices of M matrix in 
SIM, we always use Cimmino’s method. We compare our results with cgne 
which is a Krylov-type method. The method cgne is sometimes also called 
cgls. This method is scaled by M'/?, i.e., using M!/?A, M!/2b instead of 
A, b. 

As we mentioned before, Theorems 2 and 3 and the strategies (22-24) 
are based on (12). Note that the convergence analysis is not based on the 
nonascending property. Since ¢, < ¢,41 and using [8, Proposition 3.3], both 
strategies (22) and (23) satisfy (12). For the the third strategy (24) we have 


(1— Gx)" 11 — CR)? > (1 = Gaga)” 11 — C8)? 
> (1 Gog)” (1 — Ott)? 


provided that 

Cea > Ge (27) 
We do not have any mathematical proof which shows (27) holds but our 
numerical tests verifies (27) where r > 1. 


Our first tests are based on the standard head phantom from [13]. We 
report some numerical tests with an example taken from the field of tomo- 
graphic image reconstruction from projections, using the SNARKO9 software 
package [6]. The phantom is discretized into 63 x 63 pixels, and 16 projections 
(evenly distributed between 0 and 174 degrees) with 99 rays per projection 
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Figure 4: Relative error histories in SIM using small phantom with different relaxation 
strategies 


are used. The resulting matrix A has dimension 1584 x 3969, so that the 
system of equations is highly underdetermined. Figure 4 shows the error his- 
tories for SIM, using the optimal fixed relaxation parameter as well as UV), V2 
and W3 strategies with noisy data. 


Based on behavior of noise error upper bound, see Figure 3, using W3 
with r = 1 gives the smallest upper bound. This fact is confirmed by Figure 
4 (left) where 5% noise is used. But using 10% noise, Figure 4 (right), leads 
to fast semi-convergence. The reason could be the large value of w in (19) 
which is eliminated in all strategies. However, the results of V3 rule with 
r = 1.5 and W, rule are proper where 10% noise is used. 


In our second example we used the (matlab-based) package AIRtools [12] 
to produce the phantom, the matrix and the right-hand side (with and with- 
out noise). We again used 5% and 10% white Gaussian noises. The phantom 
is now discretized into 365 x 365 pixels. We take 88 projections (evenly 
distributed between 0 and 179 degrees) with 516 rays per projection. The 
resulting projection matrix A has dimension 40892 x 133225, so that again 
the system of equations is underdetermined. Figure 4 shows the relative er- 
ror histories of SIM with noisy data. As it is seen, this figure shows that the 
results of V3 rule with r = 1 are close to the results of optimal rule. 


For both noise levels and phantoms, cgne is the fastest method. However it 
also shows a distinctive semi-convergence behavior making it more dependent 
on a reliable stopping rule than SIM with our strategies. 


Acknowledgments 


We thank Tommy Elfving and Per Christian Hansen for their valuable com- 
ments. We wish to thank two anonymous referees for constructive criticism 
and helpful suggestions which improved our paper. 


Control 


ing semi-convergence phenomenon in non-stationary ... 


5% noise 


: ; ; 
—— 1, 8242 
‘ont 


aaa 
=, 
with 1 


—— with 1.5 


10% noise 


: : , : 
— i, 610 df 


Secmenll / 
=, r 
wi 1 / 


wih 21 5 
Ye with 2 ya 


63 


EY with 2 


L L L L L L L 
0 5 10 16 20 5 30 6 49 0 
iteration iteration 


Figure 5: Relative error histories in SIM using the big phantom 
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